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We study the influence of the system size on various static and dynamic properties of a super- 
cooled binary Lennard- Jones liquid via computer simulations. In this way, we demonstrate that the 
treatment of systems as small as TV = 65 particles yields relevant results for the understanding of 
bulk properties. Especially, we find that a system of TV = 130 particles behaves basically as two 
non-interacting systems of half the size. 



I. INTRODUCTION. 



The theoretical understanding of the glass transition is still far from being complete. During the last years, though, 
considerable progress has been made both from the analytical and the numerical side (see, e.g., reference [jyV. In 
this paper, we will dwell on some questions that arise within the framework of the energy landscape approach (5[ [j). 
The beauty of this viewpoint lies in the fact that the complicated nature of many-particle effects in structurally 
disordered matter can be formulated in a pictorial way by considering the topology of the high-dimensional landscape 
of the total potential energy V(x) (PEL). (Here, x = (xx,...,Xjy) denotes the set of positions of all N particles 
in the system.) Of special importance at low temperatures are the local minima of the PEL. They are extremely 
numerous, so that a statistical treatment of their properties is appropriate. From the statistics of their energies 
and vibrational characteristics, the whole thermodynamics can be derived at low-enough temperatures and constant 
volume 1. Recently, generalizations to the constant-pressure situation and non-equilibrium conditions have been 
numerically implemented [||. The critical temperature T c of mode-coupling theory serves as a good indicator 
for the temperature range where the PEL standpoint is appropriate: Below T c (the so-called landscape-dominated 
regime), it is generally accepted that the temporal evolution of a system happens through activated jumps among 
PEL minima. Between T c and 2T C (the landscape-influenced regime), properties of minima are generally deemed to 
be relevant for the thermodynamic description, whereas they are expected to be irrelevant for dynamics there. This 
has been concluded from the analysis of higher-order stationary points, which start to be populated above T c @, [§. 
In two recent publications, however, we have provided evidence that this notion should be revised: From a detailed 
analysis of relaxation dynamics, we have found that the picture of activated hopping is correct even above T c [To| . 
In any event, above 2T C , the PEL description breaks down, due to the fact that the system no longer occupies the 
well-behaved vicinity of minima. 

In our recent publications, we elucidated the implications of local PEL topology for dynamics ||, [l(], O. As conjec- 
tured by Stillinger ||, we found that the PEL is composed of groups of correlated minima, called metabasins (MB). 
Diffusional motion then turned out to consist of random jumps among metabasins, where the jump times could be 
related to the depths of metabasins. These conclusions were drawn from computational studies of fairly small model 
systems of Lennard- Jones type (see reference |10 for details about the systems). The motivation for using systems 
as small as N = 65 particles is two-fold: Firstly, much longer time scales are accessible in the simulations, and more 
sophisticated PEL analyses become possible |l(J] . Secondly, the global PEL viewpoint implies that the hypersurface of 
the potential energy is the more complex the larger the system is. On the other hand, generally, relaxation processes 
are spatially localized. This implies that PEL complexity in large systems originates mainly from a superposition of 
independently relaxing subsystems. In contrast, the local relaxation dynamics itself is governed by a non-trivial kind 
of PEL complexity which, however, is essentially identical in all the subsystems. Since we are interested in the physics 
behind local relaxation, we concentrate on very small systems. We thus avoid to consider many relaxing subsystems 
in parallel which would average out the information about a single one. 

In the present paper, we shall provide evidence that the results obtained in reference |^|, ^| for a small binary 
Lennard-Jones system of 65 particles (BMLJ65) are also relevant for the bulk behavior. Systems of N = 130, 260, and 
1000 particles will be investigated and compared to the BMLJ65. Especially, we shall demonstrate that the BMLJ130 
behaves essentially like two non-interacting BMLJ65s. 

We shall study static quantities in section || and turn to dynamic observables in section [II . Further aspects of our 



results are discussed in section [V 



II. STATIC PROPERTIES. 



Pair- correlation function. A first test for finite-size effects is to compare the distributions of interparticle distances 



2 



5 



4 



3 



2 



1 





12 3 4 

FIG. 1: Pair-correlation function, gAA(r), between A particles, for N = 1000, 130, and 65. Periodic images of the simulation 
box have been used to compute gAA(r) for distances larger than half the box width, r > Ljv/2 (N = 130 and 65). 



5a/3( r ) f° r different system sizes. Here, we restrict ourselves to the pair-correlation function of the ^4-particles, qaa{t), 
see figure [j]. Within a simulation box of width Ln, we may only calculate g a p(r) for r < Ln/2. For larger values 
of r, periodic images of the simulation box must be used. We find that gAA{f) of the BMLJ65 is identical to the 
one of the BMLJ1000 for r < L^/2. At larger distances, deviations from the bulk distribution can be seen. This is 
plausible, since the simple duplication of the simulation box can surely not reproduce all details of the long-ranged 
bulk correlations. Nevertheless, the oscillations corresponding to higher-order neighbor shells in the BMLJ1000 are 
also present in the duplicated BMLJ65. Similarly, the BMLJ130 matches the BMLJ1000 gAA{r) for r < L130/2, 
whereas deviations for larger r already seem to be negligibly small. 

Statistics of minima. It is an important question how the properties of the PEL are affected by changes in system 
size. The most prominent characteristics of a PEL minimum are its energy e and vibrational partition function 
7 1 ( 3Ar -3)/ 2 y The latter can be calculated within harmonic approximation, 




where the A„ are the eigenvalues of the hessian matrix in the minimum. Since the number of PEL minima is extremely 
large, a statistical treatment is needed. As a starting point, we analyze the mean energy of minima at temperature 
T, (e(T;N)), and their variance a 2 (T;N). For systems composed of independent subsystems, (e(T;N)/N) and 
a 2 (T;N)/N do not depend on system size. In figure ||, these quantities are shown for N = 65 and N = 130, plus 
some data points oi N — 260 and N — 1000. Concerning the mean energies, we find a good overall agreement of 
different system sizes. The maximum difference is about 1% between the BMLJ65 and the BMLJ260 at T — 00. In 
the landscape- influenced regime below T = 2T C , data for different TV show a perfect match. A similar conclusion can 
be drawn from figure ||(b), where we see a 2 (T)/N. A systematically larger value is found for the BMLJ65 at high 
temperatures, as compared to the BMLJ130. For T < 2T C , the difference is less than 20%, but more precise statements 
are prohibited by the statistical uncertainty of a 2 (T)/N below T = 0.6. Thus, small but significant finite-size effects 
can be observed in this quantity. 

We shortly comment on the deviations from the gaussian prediction at T < 0.5, as seen in figure ||(a),(b). A possible 
explanation is that our simulations below T = 0.5 were too short to sample the PEL thoroughly at the lowest energies. 
In view of the total simulation times for the BMLJ65 (T = 0.4: 50r Q , T = 0.435: 520r Q , T = 0.45: 850r Q ) this is 
strong statement. We currently study if significantly longer simulations could alter the situation in figure |^, T < 0.5. 

Total number of minima. The overall statistics of PEL minima can be described to a large degree by the number 
density of their energies, G(e). Formally, it can be written as G(e) = 8(e — ej), whereas in practice, some coarse 
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FIG. 2: (a) Mean minimum energy per particle vs. 1/T. Data of different system sizes are given, (b) Variance of minimum 
energies, also on a per-particle base, again vs. 1/T. In (a) and (b), the straight lines are the predictions from a gaussian density 
of minima. At 1/T — 0, minimizations were performed from configurations with random particle positions. 



graining of energy is introduced. In computer simulations, G(e) was found to be approximately gaussian for several 
model systems |l^, 13, 15 , 



G(e) 



V2i 



: exp aN 



-«) : 



(2) 



The lower cutoff energy e min takes into account that the ideal gaussian cannot extend to arbitrarily low energies. A 
cutoff at the high-energetic tail of G(e) is not needed, due to the small Boltzmann weight of these minima. For huge 
systems, which, to a good approximation, are composed of many independent subsystems, the gaussianity of G(e) 
is an immediate consequence of the central limit theorem. However, as discussed in reference |l6|, a large degree of 
gaussianity must already be present in the subsystems that are considered elementary. 

The modification of the gaussian by the cutoff e min is normally small, so that one finds Nq(N) = J deG(e) » e aN . 
Thus, the so-called growth parameter a describes how the total number of minima, Nq(N), evolves with system 
size. We will now calculate AT (65) and JVq(130). To this end, we apply a practical method which has recently 
been discussed in the literature |17], [l8|, [l9|, g(| ^l|. The main step is to compute the partition function Z(T) or, 
equivalently, the entropy via thermodynamic integration from a known reference state. The knowledge of Z(T) can 
be used to compute G(e) as follows pl| . Using the harmonic approximation of basin vibrations, the population of 
minimum i at temperature T is 



Z{T) 



where (3 = 1/k^T. By computing the expectation value ^ Y i 1 e /3ei J(e — ei)pi, we can then extract G(e), 

In G(e') = In (d(e - e^Y^e^) + In Z(T) - In T. 



(3) 



(4) 



In the above cited works, starting from a high-temperature (To), low-density state (N/Vq), one compressed the 
system until the required volume V\ of the supercooled liquid was reached. Subsequently, one cooled along the 
isochore down to T = T\. The procedure is depicted in figure [|(b), inset. The partition function at the state point 
(Vi, Xi) follows with the help of the relations 
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In Z(y x ,Tx) = In Z(V , T ) + f3 J dV P (V, T ) - J d0E(Vi,T). 
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FIG. 3: Determination of G(e) via thermodynamic integration, (a) Excess pressure p ox = p — NksT /V over V/N in a double- 
logarithmic plot. The straight line corresponds to the first correction to ideal-gas behavior, described by the second virial 
coefficient B2(T). (b) Temperature dependence of the mean potential energy E(T)/N = (V(x)) /N. Lines are fits of the form 
E(T) = a + bT 3 ^ 5 . Note that the data of BMLJ65 and BMLJ130 practically coincide. Inset: thermodynamic integration path 
in the V — T plane (To = 5.0, N/Vj = 1.2). (c) Number density of minimum energies, G(e), computed via equation ^ from 
simulation runs of the BMLJ65 at T = 0.4,0.435,0.45,0.466,0.48,0.5, and 0.6. (d) Number density G(e) of the BMLJ130, 
computed from T = 0.4, 0.435, 0.45, 0.5, and 0.6. In (c) and (d), data for T > 0.8 (+) do not fall onto the master curve, indicating 
the beginning breakdown of the harmonic approximation. Inset: number density of minima of BMLJ65 and BMLJ130, on a 
per-particle base. Curves coincide completely. 



In the limit Vq — ► oo, we may replace In Z(Vo, To) in equation ^ by the ideal gas term. Note that we write E(V,T) 
instead of (V) for potential energy here, in order to avoid confusion with volume. 

After measuring the pressure- volume dependence at T = 5 (see figure ||(a)), we evaluate the first integral of 
equation ^. This yields the partition function at To = 5, p\ — 1.2, 

hxZ{V u T ) =-108.3 ±0.7 (BMLJ65), (6) 
In Z(V U T ) =-218.5 ±2.5 (BMLJ130). (7) 

Considering the second integral of equation ||, we need the mean potential energy E(T) along the isochore V = V\ 
(see figure [§(b)). As is done in the literature, we use the functional form E(T) = a + bT 3 ^ 5 to parametrize our data. 
For the theoretical background of this form, see reference |22j . From the two fit parameters 

a = -361.4 ±0.3, b =171.5 ±0.3 (BMLJ65), (8) 
a = -718.8 ±0.8, 5 = 341.2 ±0.9 (BMLJ130), (9) 

one can then calculate the second integral in equation ^. 

We now use these results to calculate No, via equation ||. The expectation value in equation || is extracted from 
regular simulations. Thus, as long as the harmonic approximation holds, we are able to calculate the absolute value 
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FIG. 4: Ratios of diffusion coefficients Dn(T) for three system sizes versus 1/T. 



of G(e). For the BMLJ65 and BMLJ130 systems, G(e) is shown in figure ||(c) and (d), respectively. At temperatures 
T < 0.6, all points for G(e) fall nicely onto a master curve, whereas for T > 0.8 the normalization does not work 
anymore, sec reference |2l]]. We then fit gaussians to the data at T < 0.6, yielding the complete G(e). The number of 
minima can now be calculated from Nq = J deG(e), 



Ao = 10 22 4±0 - 8 (BMLJ65), 
No = 10 45 0±25 (BMLJ130). 



(10) 
(11) 



Hence, within error bars, we find Ao(130) = (Ao(65)) 2 , which is the trivial scaling behavior expected from combina- 
tions of non-interacting subsystems. 

We finally note that we reached the same conclusion after calculating the configurational entropy S C (T;N) as 
described in reference Q (data not shown here): S C (T; N)/N turns out to be identical for N = 65 and N = 130. 



III. DYNAMIC PROPERTIES. 



We now discuss the influence of system size on dynamics. Here, more drastic effects than in static quantities are to 
be expected: In fact, it is the most puzzling feature of the glass transition itself that a dramatic slowdown of molecular 
motion cannot be traced back to changes in static quantities easily. 

Diffusion coefficients. We start with the long-time diffusion coefficient Dn(T), defined by the Einstein relation. 
One finds that the D N (T) for N = 65,130, and 1000 differ only very little. In figure || we see -D65/-D1000 and 
Dq$/Di3q as functions of temperature. The difference between and Diooo -which we assume to be identical to 
the bulk diffusion coefficient- is twenty percent or less above T c . Since data for -Diooo are not available below T c , no 
such comparison is possible there. As reflected by D 65 /D 130 , these differences are already present between N — 65 
and N = 130. Below T c , however, the BMLJ65 seems to become slightly faster than the BMLJ130. However, since 
error bars are large for the two low-temperature data points, it is hard to judge whether this is a systematic effect 
that further increases upon cooling. In any event, in the temperature range studied, the overall variation of Dn(T) 
is more than three orders of magnitude. Regarding the small deviation of the BMLJ65 relative to the BMLJ1000, 
finite-size effects in the long-time diffusion should be jugded small. As a comparison, major changes happen when 
going to N = 40, where we find L> 40 (0.5)/L>iooo(0.5) » 0.1. 

Waiting-time distributions. As a more refined comparison of dynamics between different system sizes, we consider 
the distributions of MB lifetimes (waiting times), see references || ^3|. A detailed description of MBs and their 
properties can be found in reference JTo| . Here we briefly describe how MB lifetimes can be obtained from a given 
simulation run. Based on an equidistant time series of minima, we resolve the elementary transitions between minima 
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FIG. 5: Distributions of MB lifetimes, p(logr) at T = 0.5. As discussed in the text, the distribution pi3o(logr) should be 
reproducable from p6s(logr) by a special kind of convolution, equation [H| The corresponding function p65®65(logr) is given 
in the figure. 



by further minimizations, accompanied by temporal interval bisection. From the series of configurations thus obtained, 
we determine the time intervals of correlated back-and- forth jumps within groups of a few minima, which are identified 
with the MB lifetimes [y]. In this way, we are able to detect MB durations ranging from one MD step to 
many millions of them. At some arbitrary time of a simulation run, the probability to be in a MB of length r is 
p(r) — ^2 i TiS(r — Tj) / 'X^ r 2! where the Tj's are the lifetimes found in the run. (Since the MB residence times span 
more than six orders of magnitude, our numerical computations will involve the distributions p(\og r) = p{t)t In 10 
rather than p(r).) The temperature dependence of p(r) will be suppressed for notational convenience. 

Guided by the idea that a BMLJ130 system is basically a duplication of two independent BMLJ65s, we may 
ask whether the distribution pi3o(r) of the larger system can be reproduced by some sort of convolution of the 
distribution pe5(r) of the smaller one. (For a combined system, MB lifetimes are defined as the periods where neither 
of the subsystems relaxes.) Indeed, after a lengthy calculation one finds for the duplicated system, 



P65®65(t) = -"^ J dT'p 65 (T') J <1t"p 65 (t") ^1 - -^-jj 

T T 

This expression can be simplified and, upon using p(\og r) , it reads 
where 



P65®65(logT) = 2p 65 (log t)/(t) + 2t/(t) [I{t)t In 10 - p 65 (log r) 



(12) 



(13) 



I(t) = / 6t'p(t') and 7(r) = / &t'p{t')/t' 



In figure |(| we show p65(logr), together with pi3o(logr) at the temperature T — 0.5. The distribution resulting 
from the duplication, P65®65(l°g'?"), is also given in the figure. It agrees nicely with P130. Thus, on the refined 
level of waiting-time distributions, we find further evidence that larger systems basically behave as consisting of non- 
interacting BMLJ65-type building blocks. Essentially, p\^ is shifted to the left with respect to P65- This is no wonder, 
because time intervals where both independent systems are inert, are generally shorter than the waiting times of a 
single system. For instance, the mean waiting times obey the relation 
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FIG. 6: Apparent activation energies E app (e M B/N) derived from mean lifetimes at fixed metabasin energy e M B (equation [L4j ). 
Data for N = 65 and N = 130 are shown versus e mB /N. The interpolation of the N — 65 data has been used to compute 
-Bapp for the union of two non-interacting BMLJ65, as described in the text. The result of this calculation is given in the figure 
(solid line). Three further duplications yield -E app for N = 1040 (dashed line). 



which can be shown with the help of equation [L2|. 

Metabasin depths. As discussed above, metabasins have turned out as the relevant structures in the PEL for 
describing the slowdown of molecular motion in supercooled liquids. In a recent paper, we reported on how the 
average lifetimes (r|e MB ; T) of MBs depend on their energies e MB |fTo[ . At some fixed e M B, we found that (r|e MB ; T) is 
Arrhcnius-like below T rs 2T C , leading to the parametrization 

(r|e MB ;T) = T (e MB )e^^ . (14) 

The apparent activation energy -E app (e M B) shows a strong dependence on e MB , as soon as one drops below e MB /N w 
—4.5. This can be seen in figure || where i? ap p(e M B) versus e MB /N is depicted both for N = 65 and A" = 130. Naturally, 
it is tempting to relate E^ pp (e MB ) to PEL structure. By a detailed investigation of escape paths and the barriers that 
are overcome when escaping the MBs, we were able to reproduce i? app (e) from the local topology of MBs. Hence, we 
were able to prove that the effective depth .E app of a MB, derived from dynamics, directly corresponds to the real 
depth of that MB as given by the energy barriers around it. Moreover, escaping from deep MBs could be shown to 
involve actived jumps even above T c . What concerns the prefactor 7~o(e MB ) of equation [l4], no such understanding in 
terms of PEL structure has been possible. Fortunately, however, To(e MB ) turned out to have a weak dependence on 
MB energy. Thus, it may in good conscience be set to a constant JIM. 

Here we concentrate on the dependence of -E app (e M B) on system size. As can be seen from figure ^, the activation 
energies E app (e MB /N) of N = 65 and N — 130 are quite close. However, the N = 130 data for e MB /N < —4.5 show 
the tendency to fall slightly below that of N = 65. We shall show that this trend can be understood again in terms 
of a simple duplication of a BMLJ65. Hence, we are interested in the combination of two independent BMLJ65 
systems. Consider a MB of energy e MB = £mb + £ mb in the combined system. Then assume that its average lifetime 
can be expressed through the lifetimes of both subsystems, i.e. 

1 ii 

(15) 



(t|6mb, ei 1B ) 65(865 (t|£mb) 65 (i"|eLiB) 65 

Averaging over the population of £m B and £mb at constant e MB then yields (t\€ mb ;T) for the combined system. The 
mean lifetimes produced in this way are again Arrhenius-like below 2T C (data not shown here). Thus, data can again 
be fitted by a function of the form of equation |l4|, yielding E app (e MB ) for the duplicated BMLJ65. The result is shown 
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in figure || for T = 0.5. Again, the artificial BMLJ65 duplication reproduces the observations for the real system of 
N = 130 particles. 

Finally, we note that further duplication of the BMLJ65 leads to an interesting result: The activation energies from 
duplication nearly fall on top of each other for N > 130 and e MB /N > —4.6. The example of sixteen non-interacting 
BMLJ65s (N = 1040) is shown in figure § 

IV. DISCUSSION AND CONCLUSIONS. 

For several static and dynamic observables, we have verified the factorization property 

BMLJ130 w BMLJ65 <g> BMLJ65 . 

The BMLJ130, in turn, seems to be close to bulk behavior. Some of the results presented here have already been 
obtained in earlier work for a very similar Lennard- Jones type system [jl2). Again, the conclusion can be drawn that 
binary Lennard- Jones systems of ca. 60 particles are a very good compromise between the desired smallness needed for 
our PEL investigations and the required absence of large finite-size related artifacts. In this connection, the simulation 
results of Yamamoto and Kim on the standard binary soft sphere mixture are of interest p4[ . Comparing systems 
of N = 108 and 1000 particles above T c , the authors found the small system to be up to an order of magnitude 
slower than the large one. These findings suggest a fundamental difference between the Lennard-Jones and soft- 
sphere types of relaxation dynamics. Evidently, soft spheres exibit a larger length scale of cooperative motion than 
do Lennard-Jones systems. A possible explanation can be found in reference [Q. 

It is known from the study of cooperative length scales that they increase with decreasing temperature [^5|, ^6], ^7) . 
Thus, at some lower temperatures one might expect that 65 particles are no longer enough and finite-size effects 
become visible. For a similar Lennard-Jones system it has been shown that finite-size effects are reflected by the fact 
that the bottom of the PEL is frequently probed jl^] . In the present case still longer simulations at lower temperatures 
have to be performed to check whether also for N = 65 the PEL bottom can be reached. Then the interesting question 
arises whether or not differences to the N = 130 system become visible. 

The results presented in this paper suggest that the essential physics of the supercooled BMLJ is already contained 
in the system of N = 65 particles. For the temperatures under investigation here, extended structures of collectively 
moving particles ('strings') have been reported in large systems [^8|. It is important to know whether these structures 
are still present in the small systems considered here. Work along this line is in progress. 
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